
#include "PyrolysisSolidADMaterial.h"

registerMooseObject("TrilobitaApp", PyrolysisSolidADMaterial);

InputParameters
PyrolysisSolidADMaterial::validParams()
{
    auto params = ADMaterial::validParams();
   
    params.addRequiredCoupledVar("rho_A", "Coupled density for pahseA");
    params.addRequiredCoupledVar("rho_B", "Coupled density for pahseB");
    params.addRequiredCoupledVar("rho_C", "Coupled density for pahseC");
    params.addParam<Real>("volume_fraction_A", 0, "volume fraction of phaseA");
    params.addParam<Real>("volume_fraction_B", 0, "volume fraction of phaseB");
    params.addParam<Real>("volume_fraction_C", 0, "volume fraction of phaseC");
    params.addRequiredCoupledVar("temperature", "Coupled temperature");
    params.addRequiredParam<Real>("rhov", "density of original material");
    params.addRequiredParam<Real>("rhoc", "density of fully charred material");
    params.addRequiredParam<std::vector<Real> >("kv_list", "The vector of kv values of origin material for building the piecewise function");
    params.addRequiredParam<std::vector<Real> >("kvT_list", "The vector of temperature values for building the piecewise function to compute kv ");
    params.addRequiredParam<std::vector<Real> >("kc_list", "The vector of kc values for building the piecewise function to compute kc ");
    params.addRequiredParam<std::vector<Real> >("kcT_list", "The vector of temperature values for building the piecewise function to compute kc ");
    params.addRequiredParam<std::vector<Real> >("cpv_list", "The vector of cpv values for building the piecewise function to compute thermal conductivity of origin material ");
    params.addRequiredParam<std::vector<Real> >("cpvT_list", "The vector of cpv values for building the piecewise function to compute thermal conductivity of origin material ");
    params.addRequiredParam<std::vector<Real> >("cpc_list", "The vector of cpc values for building the piecewise function to compute cpc ");
    params.addRequiredParam<std::vector<Real> >("cpcT_list", "The vector of temperature values for building the piecewise function to compute cpc");
    params.addRequiredParam<Real>("pyrolysis_heat", "pyrolysis heat");
    params.addRequiredParam<Real>("porosity_v", "porosity of original material");
    params.addRequiredParam<Real>("porosity_c", "porosity of fully charred material");
    params.addRequiredParam<Real>("permeability_v", "permeability of original material");
    params.addRequiredParam<Real>("permeability_c", "permeability of fully charred material");
    params.addRequiredParam<Real>("pyrolysis_heat", "pyrolysis heat");
    params.addRequiredParam<Real>("surface_emissivity", "emissivity");

    return params;
}

PyrolysisSolidADMaterial::PyrolysisSolidADMaterial(const InputParameters& parameters)
    : ADMaterial(parameters),
    _rho_A_value(adCoupledValue("rho_A")),
    _rho_B_value(adCoupledValue("rho_B")),
    _rho_C_value(adCoupledValue("rho_C")),
    _rho_A_dt(adCoupledDot("rho_A")),
    _rho_B_dt(adCoupledDot("rho_B")),
    _rho_C_dt(adCoupledDot("rho_C")),
    _temperature(coupledValue("temperature")),
    _rhos(declareADProperty<Real>("solid_density")),
    _drhos_dt(declareADProperty<Real>("drhos_dt")),
    _rhov(declareADProperty<Real>("rhov")),
    _rhoc(declareADProperty<Real>("rhoc")),
    _char_percent(declareADProperty<Real>("char_percent")),
    _thermal_conductivity_v(declareADProperty<Real>("solid_thermal_conductivity_v")),
    _thermal_conductivity_c(declareADProperty<Real>("solid_thermal_conductivity_c")),
    _thermal_conductivity_s(declareADProperty<Real>("solid_thermal_conductivity")),
    _specific_heat_v(declareADProperty<Real>("solid_specific_heat_v")),
    _specific_heat_c(declareADProperty<Real>("solid_specific_heat_c")),
    _specific_heat_s(declareADProperty<Real>("solid_specific_heat")),
    _porosity_v(declareADProperty<Real>("porosity_v")),
    _porosity_c(declareADProperty<Real>("porosity_c")),
    _porosity(declareADProperty<Real>("porosity")),
    _permeability_v(declareADProperty<Real>("permeability_v")),
    _permeability_c(declareADProperty<Real>("permeability_c")),
    _permeability(declareADProperty<Real>("permeability")),
    _pyrolysis_heat(declareADProperty<Real>("pyrolysis_heat")),
    _surface_emissivity(declareProperty<Real>("surface_emissivity"))
    
{
    _volume_fraction_A_value = getParam<Real>("volume_fraction_A");
    _volume_fraction_B_value = getParam<Real>("volume_fraction_B");
    _volume_fraction_C_value = getParam<Real>("volume_fraction_C");
    _rhov_value = getParam<Real>("rhov");
    _rhoc_value = getParam<Real>("rhoc");
    _kvT_list = getParam<std::vector<Real>>("kvT_list");
    _kv_list = getParam<std::vector<Real>>("kv_list");
    _kcT_list = getParam<std::vector<Real>>("kcT_list");
    _kc_list = getParam<std::vector<Real>>("kc_list");
    _cpvT_list = getParam<std::vector<Real>>("cpvT_list");
    _cpv_list = getParam<std::vector<Real>>("cpv_list");
    _cpcT_list = getParam<std::vector<Real>>("cpcT_list");
    _cpc_list = getParam<std::vector<Real>>("cpc_list");
    _func_kv_T = LinearInterpolation(_kvT_list, _kv_list);
    _func_kc_T = LinearInterpolation(_kcT_list, _kc_list);
    _func_cpv_T = LinearInterpolation(_cpvT_list, _cpv_list);
    _func_cpc_T = LinearInterpolation(_cpcT_list, _cpc_list);
    _pyrolysis_heat_value = getParam<Real>("pyrolysis_heat");
    _porosity_v_value = getParam<Real>("porosity_v");
    _porosity_c_value = getParam<Real>("porosity_c");
    _permeability_v_value = getParam<Real>("permeability_v");
    _permeability_c_value = getParam<Real>("permeability_c");
    _pyrolysis_heat_value = getParam<Real>("pyrolysis_heat");
    _surface_emissivity_value = getParam<Real>("surface_emissivity");
}

void
PyrolysisSolidADMaterial::initQpStatefulProperties()
{
    ADReal rhoscom = _rho_A_value[_qp] * _volume_fraction_A_value + _rho_B_value[_qp] * _volume_fraction_B_value + _rho_C_value[_qp] * _volume_fraction_C_value;
    if ((rhoscom - _rhoc[_qp]) > 0)
    {
        _rhos[_qp] = rhoscom;
    }
    else
    {
        _rhos[_qp] = _rhoc[_qp];
    }
    if ((rhoscom - _rhov[_qp]) > 0)
    {
        _rhos[_qp] = _rhov[_qp];
    }
    else
    {
        _rhos[_qp] = rhoscom;
    }
    _drhos_dt[_qp] = _rho_A_dt[_qp] * _volume_fraction_A_value + _rho_B_dt[_qp] * _volume_fraction_B_value + _rho_C_dt[_qp] * _volume_fraction_C_value;
    _rhov[_qp] = _rhov_value;
    _rhoc[_qp] = _rhoc_value;
    _thermal_conductivity_v[_qp] = _func_kv_T.sample(_temperature[_qp]);
    _thermal_conductivity_c[_qp] = _func_kc_T.sample(_temperature[_qp]);
    _specific_heat_v[_qp] = _func_cpv_T.sample(_temperature[_qp]);
    _specific_heat_c[_qp] = _func_cpc_T.sample(_temperature[_qp]);
    _char_percent[_qp] = (_rhov_value - _rhos[_qp]) / (_rhov_value - _rhoc_value + 1e-8);
    _thermal_conductivity_s[_qp] = _thermal_conductivity_v[_qp] * (1 - _char_percent[_qp]) + _thermal_conductivity_c[_qp] * _char_percent[_qp];;
    _specific_heat_s[_qp] = _specific_heat_v[_qp] * (1 - _char_percent[_qp]) + _specific_heat_c[_qp] * _char_percent[_qp];
    _porosity_v[_qp] = _porosity_v_value;
    _porosity_c[_qp] = _porosity_c_value;
    _permeability_v[_qp] = _permeability_v_value;
    _permeability_c[_qp] = _permeability_c_value;
    _porosity[_qp] = _porosity_v[_qp] * (1 - _char_percent[_qp]) + _porosity_c[_qp] * _char_percent[_qp];
    _permeability[_qp] = _permeability_v[_qp] * (1 - _char_percent[_qp]) + _permeability_c[_qp] * _char_percent[_qp];
    _pyrolysis_heat[_qp] = _pyrolysis_heat_value;
    _surface_emissivity[_qp] = _surface_emissivity_value;
}


void
PyrolysisSolidADMaterial::computeQpProperties()
{
    ADReal rhoscom = _rho_A_value[_qp] * _volume_fraction_A_value + _rho_B_value[_qp] * _volume_fraction_B_value + _rho_C_value[_qp] * _volume_fraction_C_value;
    if ((rhoscom - _rhoc[_qp]) > 0)
    {
        _rhos[_qp] = rhoscom;
        _drhos_dt[_qp] = _rho_A_dt[_qp] * _volume_fraction_A_value + _rho_B_dt[_qp] * _volume_fraction_B_value + _rho_C_dt[_qp] * _volume_fraction_C_value;
    }
    else
    {
        _rhos[_qp] = _rhoc[_qp];
        _drhos_dt[_qp] = 0;
    }
    
    _rhov[_qp] = _rhov_value;
    _rhoc[_qp] = _rhoc_value;
    _thermal_conductivity_v[_qp] = _func_kv_T.sample(_temperature[_qp]);
    _thermal_conductivity_c[_qp] = _func_kc_T.sample(_temperature[_qp]);
    _specific_heat_v[_qp] = _func_cpv_T.sample(_temperature[_qp]);
    _specific_heat_c[_qp] = _func_cpc_T.sample(_temperature[_qp]);
    _char_percent[_qp] = (_rhov_value - _rhos[_qp]) / (_rhov_value - _rhoc_value + 1e-8);
    _thermal_conductivity_s[_qp] = _thermal_conductivity_v[_qp] * (1 - _char_percent[_qp]) + _thermal_conductivity_c[_qp] * _char_percent[_qp];;
    _specific_heat_s[_qp] = _specific_heat_v[_qp] * (1 - _char_percent[_qp]) + _specific_heat_c[_qp] * _char_percent[_qp];
    _porosity_v[_qp] = _porosity_v_value;
    _porosity_c[_qp] = _porosity_c_value;
    _permeability_v[_qp] = _permeability_v_value;
    _permeability_c[_qp] = _permeability_c_value;
    _porosity[_qp] = _porosity_v[_qp] * (1 - _char_percent[_qp]) + _porosity_c[_qp] * _char_percent[_qp];
    _permeability[_qp] = _permeability_v[_qp] * (1 - _char_percent[_qp]) + _permeability_c[_qp] * _char_percent[_qp];
    _pyrolysis_heat[_qp] = _pyrolysis_heat_value;
    _surface_emissivity[_qp] = _surface_emissivity_value;
}